1 Setup R

rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)

library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(deSolve)
library(rootSolve)

set.seed(13)

run_sims <- FALSE

2 Introduction

The aim is to reproduce the model based results / figures.

Here is the original publication: https://www.nature.com/articles/s41467-017-00912-x

And here is the supplementary information, including the table of parameter values: https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-017-00912-x/MediaObjects/41467_2017_912_MOESM1_ESM.pdf

Note that below I go with the notation from the table S1 and not the main text. However there is only one small difference between these: In the main text the notation is simplified to y_S_PB (which is y_SR_PB in Table S1) and y_S_SB (which is y_SO_SB in Table S1).

The published model is modified:

  • Random noise in substrate concentrations can be included.
  • Temporal changes in oxygen diffusitivity can be included.
  • Immigration can be included.
  • A minimum abundance of each species can be set.

3 Methods

3.1 Model dimensions

  • Time: hours
  • Volume: litres
  • Substrate quantity: micromoles
  • Organism quantity: cells

3.2 Dynamic model

source("functions/bushplus_dynamic_model.r")

3.3 Defaults

3.3.1 Default parameter values

default_parameter_values <- c(
  
  ## maximum specific growth rates
  g_max_CB = 0.05,
  g_max_PB = 0.07,
  g_max_SB = 0.1,
  
  ## half-saturation constants
  k_PB_SR = 10,
  k_SB_SO = 5,
  k_CB_P = 0.2,
  k_PB_P = 0.5,
  k_SB_P = 0.5,
  
  ## half-inhibition constant
  h_SR_CB = 300,
  h_O_PB = 100,
  h_O_SB = 100,
  
  ## yield (cells per micromole)
  Y_SO_SB = 3.33e7,  
  Y_SR_PB = 1.25e7, ## Y_S_PB in the main text
  y_P_CB = 1.67e8,
  y_P_PB = 1.67e8,
  y_P_SB = 1.67e8,
  
  ## oxygen production per microbial cell
  P_CB = 6e-9,
  
  ## mortality rate
  m_CB = 0.02,
  m_PB = 0.028,
  m_SB = 0.04,
  
  ## substrate diffusivity
  a_S = 0.001,
  a_O = 8e-4,
  a_P = 0.01,
  
  ## background substrate concentration
  back_SR = 300,
  back_SO = 300,
  back_O = 300,
  back_P = 9.5,
  
  ## immigration rates (cells per time)
  i_CB <- 0,
  i_PB <- 0,
  i_SB <- 0,
  
  ## oxidisation rate of reduced sulphur
  c = 4e-5
  )

3.3.2 Default event

3.3.2.1 Default event definition

Used to add noise and to ensure minimum abundance. Noise is added to the resource concentrations, and not to the organism abundances. Minimum organism abundance is set. I.e. organismal abundance cannot go below some very low level.

source("functions/default_event_definition.r")

3.3.2.2 Default event schedule

This is when an event happens

default_event_interval <- 10

3.3.2.3 Default noise event options

Default to no noise:

default_noise_sigma <- 0

3.3.2.4 Default minimum abundance event

Default to following minimums

default_minimum_abundances <- rep(1e2, 3)
names(default_minimum_abundances) <- c("CB", "PB", "SB") 

3.3.3 Default simulation duration

default_sim_duration <- 2000
default_sim_sample_interval <- 1

3.3.4 Default oxygen diffusivity (a) forcing

Here we specify no forcing of a_0 (used unless set elsewhere).

default_log10a_series <- c(log10(default_parameter_values["a_O"]),
                           log10(default_parameter_values["a_O"]))

3.3.5 Default initial state

Initial conditions favouring anoxic

default_initial_state <- c(N_CB = 5e1,
           N_PB = 1e7,
           N_SB = 1e7,
           SO = 300,
           SR = 300,
           O = 1e1,
           P = 1e1)

3.4 Simulation function

source("functions/run_simulation.r")

3.5 Helper functions

source("functions/plot_dynamics.r")
source("functions/ss_expt_a_N.r")
source("functions/get_nonlinearity.r")
source("functions/get_hysteresis_range.r")

3.6 Code test

# dynamic_model = bushplus_dynamic_model
# event_definition = default_event_definition
# parameter_values = default_parameter_values
# event_interval = default_event_interval
# noise_sigma = default_noise_sigma
# minimum_abundances = default_minimum_abundances
# sim_duration = default_sim_duration
# sim_sample_interval = default_sim_sample_interval
# log10a_forcing = default_log10a_forcing
# initial_state = default_initial_state
# solver_method = "radau"
# 
# 
# times <- seq(0,
#              sim_duration,
#              by = sim_sample_interval)
# 
# event_times <- seq(min(times),
#                    max(times),
#                    by = event_interval)  
# 
# log10a_forcing <- approxfun(x = log10a_forcing[,1],
#                             y = log10a_forcing[,2],
#                             method = "linear",
#                             rule = 2)
# assign("log10a_forcing", log10a_forcing, envir = .GlobalEnv)
# 
# out <- as.data.frame(ode(y = initial_state,
#                          times = times,
#                          func = dynamic_model,
#                          parms = parameter_values,
#                          method = solver_method,
#                          events = list(func=event_definition,
#                                        time=event_times)))
# 
# sim_res <- run_simulation(dynamic_model = bushplus_dynamic_model,
#                            event_definition = default_event_definition,
#                            parameter_values = default_parameter_values,
#                            event_interval = default_event_interval,
#                            noise_sigma = default_noise_sigma,
#                            minimum_abundances = default_minimum_abundances,
#                            sim_duration = default_sim_duration,
#                            sim_sample_interval = default_sim_sample_interval,
#                            log10a_forcing = default_log10a_forcing,
#                            initial_state = default_initial_state,
#                            solver_method = "radau")

#sim_res <- run_simulation()

#plot_dynamics(sim_res)

4 Results: Reproducing previous results

4.1 Initial conditions favouring anoxic, no noise, no forcing

Figure 2a, b in Bush et al.

init_state_favouring_anoxic <- c(N_CB = 5e1,
           N_PB = 1e7,
           N_SB = 1e7,
           SO = 300,
           SR = 300,
           O = 1e1,
           P = 1e1)

sim_res <- run_simulation(initial_state = init_state_favouring_anoxic)
plot_dynamics(simulation_result = sim_res)

#ggsave('bacteria_no_noise.png')

Population dynamics seem as in the paper. Resource dynamics fluctuate more than in the paper.

4.2 Initial conditions favouring anoxic, with noise, no forcing

Figure 2a, b in Bush et al, but with noise

## add some noise
noise_sigma_0.01 <- 0.01
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic,
                          noise_sigma = noise_sigma_0.01)
plot_dynamics(sim_res)

#ggsave("high_SB_PB.png")

4.3 Initial conditions favouring oxic, no noise, no forcing

Figure 2c, d in Bush et al.

init_state_favouring_oxic <- c(N_CB = 1e7,
           N_PB = 1e5,
           N_SB = 1e5,
           SO = 500,
           SR = 50,
           O = 30,
           P = 10)



sim_res <- run_simulation(initial_state = init_state_favouring_oxic)

plot_dynamics(sim_res)

##ggsave("high_CB.png")

Population and resource dynamics seen as those in the paper.

4.4 Initial conditions favouring oxic, with noise, no forcing

Figure 2c, d in Bush et al. but with noise.

sim_res <- run_simulation(initial_state = init_state_favouring_oxic,
                          noise_sigma = noise_sigma_0.01)

plot_dynamics(sim_res)

##ggsave("high_CB.png")

5 Results: Composition/diversity and nature of response

ssfind_minimum_abundances <- rep(0, 3)
names(ssfind_minimum_abundances) <- c("CB", "PB", "SB") 
ssfind_simulation_duration <- 10000
ssfind_simulation_sampling_interval <- ssfind_simulation_duration
ssfind_event_interval <- ssfind_simulation_duration
ssfind_parameters <- default_parameter_values
# ssfind_init_state <- c(N_CB = NA, 
#                        N_PB = NA,
#                        N_SB = NA,
#                        SO = 500,
#                        SR = 50,
#                        O = 30,
#                        P = 10)
##ssfind_init_state <- default_initial_state
ssfind_init_state <- init_state_favouring_anoxic
ssfind_init_state <- c(N_CB = NA, #this one is varied
           N_PB = 1e8,
           N_SB = 1e8,
           SO = 200,
           SR = 200,
           O = 10,
           P = 9.5)
grid_num_a <- 100
a_Os <- 10^seq(-4.5, -1.5, length=grid_num_a)
grid_num_N <- 2

5.0.1 Each composition

5.0.1.1 No organisms

initial_CBs <- 0
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)


ss_no_biology <- ss_by_a_N(expt)
saveRDS(ss_no_biology, "data/ss_no_biology.RDS")
ss_no_biology <- readRDS("data/ss_no_biology.RDS")
ss_no_biology %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_no_biology %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.2 CB

initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)


ss_only_CB <- ss_by_a_N(expt)
saveRDS(ss_only_CB, "data/ss_only_CB.RDS")
ss_only_CB <- readRDS("data/ss_only_CB.RDS")
tempx1 <- ss_only_CB %>%
  select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.))
tempx1 %>%
  ggplot(aes(x = a, y = log10(Quantity), col = Species)) +
  geom_path(size=0.8) +
  ylab(expression(Log[10](Abundance))) +
  xlab(expression(Log[10](Oxygen~diffusivity))) +
  theme_light() +
  ylim(0,10) +
  geom_path(data = filter(tempx1, Species == "PB"),
            aes(y=0.2), size = 0.8) +
  theme(plot.title=element_text(size=10, hjust=0.5)) +
  labs(title="Only CB") 

ggsave("figures/CB.png", width = 3, height = 2)
ss_only_CB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.3 SB

initial_CBs <- 0  # this one is varied
initial_PBs <- 0
initial_SBs <- 10^seq(0, 8, length=grid_num_N)
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)


ss_only_SB <- ss_by_a_N(expt)

saveRDS(ss_only_SB, "data/ss_only_SB.RDS")
ss_only_SB <- readRDS("data/ss_only_SB.RDS")
ss_only_SB %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_only_SB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.4 PB

initial_CBs <- 0  # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)


ss_only_PB <- ss_by_a_N(expt)

saveRDS(ss_only_PB, "data/ss_only_PB.RDS")
ss_only_PB <- readRDS("data/ss_only_PB.RDS")
tempx1 <- ss_only_PB %>%
  select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.))
tempx1 %>%
  ggplot(aes(x = a, y = log10(Quantity), col = Species)) +
  geom_path(size = 0.8) +
  theme_light() +
  ylim(0,10) +
  ylab(expression(Log[10](Abundance))) +
  xlab(expression(Log[10](Oxygen~diffusivity))) +
  geom_path(data = filter(tempx1, Species == "CB"),
            aes(y=0.2), size = 0.8)+
  geom_path(data = filter(tempx1, Species == "PB"),
            size = 0.8) +
  theme(plot.title=element_text(size=10, hjust=0.5)) +
  labs(title="Only PB") 

ggsave("figures/PB.png", width = 3, height = 2)
ss_only_PB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.5 SB and PB

initial_CBs <- 0  # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)
expt$N_SB <- expt$N_PB

ss_only_SB_PB <- ss_by_a_N(expt)
saveRDS(ss_only_SB_PB, "data/ss_only_SB_PB.RDS")
ss_only_SB_PB <- readRDS("data/ss_only_SB_PB.RDS")
ss_only_SB_PB %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_only_SB_PB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.6 SB and CB

initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)
expt$N_SB <- expt$N_CB

ss_only_SB_CB <- ss_by_a_N(expt)
saveRDS(ss_only_SB_CB, "data/ss_only_SB_CB.RDS")
ss_only_SB_CB <- readRDS("data/ss_only_SB_CB.RDS")
ss_only_SB_CB %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_only_SB_CB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:5) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.7 CB and PB

initial_CBs <- 0  # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)
expt$N_CB <- expt$N_PB

ss_only_CB_PB <- ss_by_a_N(expt)
saveRDS(ss_only_CB_PB, "data/ss_only_CB_PB.RDS")
ss_only_CB_PB <- readRDS("data/ss_only_CB_PB.RDS")
ss_only_CB_PB %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_only_SB_PB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.1.8 CB, SB, and PB

initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 1e5
initial_SBs <- 1e5
expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    a_O = a_Os)

ss_CB_SB_PB <- ss_by_a_N(expt)
saveRDS(ss_CB_SB_PB, "data/ss_CB_SB_PB.RDS")
ss_CB_SB_PB <- readRDS("data/ss_CB_SB_PB.RDS")
ss_CB_SB_PB %>%
  select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species),
            size=0.8) +
  ylab(expression(Log[10](Abundance))) +
  xlab(expression(Log[10](Oxygen~diffusivity))) +
  theme_light() +
  ylim(0,10) +
  theme(plot.title=element_text(size=10, hjust=0.5)) +
  labs(title="CB, SB, & PB") 

ggsave("figures/CB_SB_PB.png", width = 3, height = 2)
ss_CB_SB_PB %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

5.0.2 All response measures

all_ss <- bind_cols(readRDS("data/ss_no_biology.RDS"),
                    composition = rep("No organisms", grid_num_a*grid_num_N),
                    direction = c(rep("up", grid_num_a),
                                  rep("down", grid_num_a)))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_only_CB.RDS"),
                              composition = rep("CB", grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_only_SB.RDS"),
                              composition = rep("SB", grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_only_PB.RDS"),
                              composition = rep("PB", grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_only_SB_PB.RDS"),
                              composition = rep("SB-PB", grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_only_SB_CB.RDS"),
                              composition = rep("SB-CB", grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_only_CB_PB.RDS"),
                              composition = rep("CB-PB", grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
                    bind_cols(readRDS("data/ss_CB_SB_PB.RDS"),
                              composition = rep("CB-SB-PB",
                                                grid_num_a*grid_num_N),
                              direction = c(rep("up", grid_num_a),
                                            rep("down", grid_num_a))))


x_long <- all_ss %>%
  gather(key = Species, value = Quantity, 2:8)
x_wide <- x_long %>%
  select(-initial_N_CB) %>%
  spread(key = direction, value=Quantity, drop=T)

resp_summary <- x_wide %>%
  group_by(composition, Species) %>%
  summarise(hyst_1 = mean(abs(log10(up) - log10(down))),
            hyst_2 = get_hysteresis_range(up, down, a),
            nl_up = get_nonlinearity(a, log10(up)),
            nl_down = get_nonlinearity(a, log10(down))) %>%
  ungroup() %>%
  gather(key = resp_measure, value = value, 3:6) %>%
  mutate(composition = fct_relevel(composition,
                                   "No organisms",
                                   "CB",
                                   "SB",
                                   "PB",
                                   "CB-PB",
                                   "SB-CB",
                                   "SB-PB",
                                   "CB-SB-PB"),
         num_func_groups = case_when(composition =="No organisms" ~ 0,
                                     composition == "CB" ~ 1,
                                     composition == "SB" ~ 1,
                                     composition == "PB" ~ 1,
                                     composition == "CB-PB" ~ 2,
                                     composition == "SB-CB" ~ 2,
                                     composition == "SB-PB" ~ 2,
                                     composition == "CB-SB-PB" ~ 3
                                     ),
         var_type = ifelse(str_sub(Species, 1, 1)=="N", "Organism", "Substrate")
         )
saveRDS(resp_summary, "data/resp_summary.RDS")
resp_summary <- readRDS("data/resp_summary.RDS")
resp_summary %>%
  ggplot() +
  geom_point(aes(x = composition, y = value)) +
  facet_grid(Species ~ resp_measure, scales = "free") +
  coord_flip()

5.0.3 Nonlinearity

5.0.3.1 All species nonlinearity

resp_summ1 <- resp_summary %>%
  filter(resp_measure %in% c("nl_down", "nl_up")) %>%
  group_by(composition, Species, num_func_groups, var_type) %>%
  summarise(ave_nl = mean(value))

resp_summ1 %>%
  ggplot() +
  geom_point(aes(composition, ave_nl)) +
  facet_wrap( ~ Species)

5.0.3.2 Substrate nonlinearity

focus_on <- "Substrate"
ave_nl_substrate <- resp_summ1 %>%
  filter(var_type == focus_on) %>%
  group_by(composition) %>%
  summarise(ave_ave_nl = mean(ave_nl))
ave_nl_substrate %>%
  ggplot() +
  geom_point(aes(composition, ave_ave_nl)) +
  ggtitle(paste("Nonlinearity in", focus_on, "response."))

5.0.3.3 Organism nonlinearity

focus_on <- "Organism"
ave_nl_organism <- resp_summ1 %>%
  filter(var_type == focus_on) %>%
  group_by(composition) %>%
  summarise(ave_ave_nl = sum(ave_nl)/unique(num_func_groups))
ave_nl_organism %>%
  ggplot() +
  geom_point(aes(composition, ave_ave_nl)) +
  ggtitle(paste("Nonlinearity in", focus_on, "response."))
## Warning: Removed 1 rows containing missing values (geom_point).

5.0.3.4 System nonlinearity

focus_on <- "System"
ave_nl_system <- resp_summ1 %>%
  group_by(composition) %>%
  summarise(ave_ave_nl = sum(ave_nl) / (unique(num_func_groups)+4))
ave_nl_system %>%
  ggplot() +
  geom_point(aes(composition, ave_ave_nl)) +
  ggtitle(paste("Nonlinearity in", focus_on, "response."))

ggsave("figures/system_nonlinearity.png", width = 3, height = 2)

5.0.4 Hysteresis

5.0.4.1 Substrate hysteresis

focus_on <- "Substrate"
resp_summ3 <- resp_summary %>%
  filter(resp_measure %in% c("hyst_1", "hyst_2"),
         var_type == focus_on) %>%
  group_by(composition, num_func_groups, resp_measure, var_type) %>%
  summarise(ave_hyst = mean(value))

resp_summ3 %>%
  ggplot() +
  geom_point(aes(x = composition, y = ave_hyst)) +
  facet_wrap( ~ resp_measure, scales="free") +
  ggtitle(paste("Hysteresis in", focus_on, "response."))

5.0.4.2 Organism hysteresis

focus_on <- "Organism"
resp_summ3 <- resp_summary %>%
  filter(resp_measure %in% c("hyst_1", "hyst_2"),
         var_type == focus_on) %>%
  group_by(composition, num_func_groups, resp_measure, var_type) %>%
  summarise(ave_hyst = sum(value) / unique(num_func_groups))

resp_summ3 %>%
  ggplot() +
  geom_point(aes(x = composition, y = ave_hyst)) +
  facet_wrap( ~ resp_measure, scales="free") +
  ggtitle(paste("Hysteresis in", focus_on, "response."))
## Warning: Removed 2 rows containing missing values (geom_point).

5.0.4.3 System hysteresis (separate)

resp_summ4 <- resp_summary %>%
  filter(resp_measure %in% c("hyst_1", "hyst_2")) %>%
  group_by(resp_measure, composition) %>%
  summarise(mean_hyst = sum(value) / (unique(num_func_groups)+4))

resp_summ4 %>%
  ggplot() +
  geom_point(aes(x = composition, y = mean_hyst)) +
  facet_wrap( ~ resp_measure) +
  ggtitle("Hysteresis in system response.")

5.0.4.4 System hysteresis (combined)

resp_summ5 <- resp_summ4 %>%
  group_by(composition) %>%
  summarise(mean_mean_hyst = mean(mean_hyst))

resp_summ5 %>%
  ggplot() +
  geom_point(aes(x = composition, y = mean_mean_hyst)) +
  xlab("Functional groups present") +
  ylab("Amount of hysteresis") +
  ylim(0,2) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggsave("figures/composition_system_hysteresis.png", width = 2.5, height = 2.5)

6 Results: Varying strength of inhibition

6.1 Introduction

Three parameters set the strength of inhibition (half-inhibition constants):

  • h_SR_CB = 300,
  • h_O_PB = 100,
  • h_O_SB = 100,

These work by the inhibition function:

\(\frac{1}{1 + \frac{X}{h_X}}\)

inhibition
## function (X, h_X) 
## 1/(1 + (X/h_X))

Plot some examples:

default_parameter_values["h_SR_CB"]
## h_SR_CB 
##     300
SR <- 10^seq(1.4, 2.5, length = 100)
ggplot() +
  geom_line(aes(x = SR, y = inhibition(SR, default_parameter_values["h_SR_CB"]))) +
  ylab("Multiplicative effect\non CB growth rate")

default_parameter_values["h_O_PB"]
## h_O_PB 
##    100
O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
  geom_line(aes(x = O, y = inhibition(O, default_parameter_values["h_O_PB"]))) +
  ylab("Multiplicative effect\non PB growth rate")

default_parameter_values["h_O_SB"]
## h_O_SB 
##    100
O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
  geom_line(aes(x = O, y = inhibition(O, default_parameter_values["h_O_SB"]))) +
  ylab("Multiplicative effect\non SB growth rate")

proportion_change <- 2
default_parameter_values["h_O_SB"]
## h_O_SB 
##    100
O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
  geom_line(aes(x = O, y = inhibition(O, proportion_change*default_parameter_values["h_O_SB"]))) +
  ylab("Multiplicative effect\non SB growth rate")

6.2 Experiment

Increasing the inhibition half saturation constants weakens the inhibition.

ssfind_minimum_abundances <- default_minimum_abundances
ssfind_simulation_duration <- 10000
ssfind_simulation_sampling_interval <- ssfind_simulation_duration
ssfind_event_interval <- ssfind_simulation_duration
ssfind_parameters <- default_parameter_values
ssfind_init_state <- c(N_CB = NA, #this one is varied
           N_PB = NA,
           N_SB = NA,
           SO = 200,
           SR = 200,
           O = 10,
           P = 9.5)

grid_num_mult <- 20
grid_num_N <- 2
grid_num_a <- 100

a_Os <- 10^seq(-4.5, -1.5, length=grid_num_a)
inhib_const_mults <- 2^seq(-1, 0.4, length=grid_num_mult)


initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 1e5
initial_SBs <- 1e5

expt <- expand.grid(N_CB = initial_CBs,
                    N_PB = initial_PBs,
                    N_SB = initial_SBs,
                    inhib_const_mult = inhib_const_mults,
                    a_O = a_Os)
ss_h1 <- ss_by_a_h_N(expt)

saveRDS(ss_h1, "data/ss_h1.RDS")
ss_h1 <- readRDS("data/ss_h1.RDS")

6.3 Results

6.3.1 Close to default inhibition

inhib_const_mults
##  [1] 0.5000000 0.5262004 0.5537737 0.5827919 0.6133306 0.6454696 0.6792927
##  [8] 0.7148882 0.7523489 0.7917725 0.8332620 0.8769256 0.9228772 0.9712366
## [15] 1.0221302 1.0756906 1.1320576 1.1913783 1.2538074 1.3195079
focal_inhibition_multiplier <- inhib_const_mults[14]

ss_h1 %>%
  filter(inhib_const_mult == focal_inhibition_multiplier) %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_h1 %>%
  filter(inhib_const_mult == focal_inhibition_multiplier) %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

6.3.2 Double strength inhibition

## YES! 0.5 does mean twice as strong inhibition
focal_inhibition_multiplier <- 0.5

ss_h1 %>%
  filter(inhib_const_mult == focal_inhibition_multiplier) %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_h1 %>%
  filter(inhib_const_mult == focal_inhibition_multiplier) %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

6.3.3 Weaker inhibition

## YES! 2 does mean half as strong inhibition
focal_inhibition_multiplier <- 2^0.4

ss_h1 %>%
  filter(inhib_const_mult == focal_inhibition_multiplier) %>%
  select(a, N_CB, N_SB, N_PB) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

ss_h1 %>%
  filter(inhib_const_mult == focal_inhibition_multiplier) %>%
  select(a, SO, SR, O, P) %>%
  gather(key = Species, value=Quantity, 2:ncol(.)) %>%
  ggplot() +
  geom_path(aes(x = a, y = log10(Quantity), col = Species))

6.3.4 Nonlinearity

x_long <- ss_h1 %>%
  mutate(direction = ifelse(initial_N_CB==1, "up", "down")) %>%
  gather(key = Species, value = Quantity, 3:9) 
x_wide <- x_long %>%
  select(-initial_N_CB) %>%
  spread(key = direction, value=Quantity, drop=T)

resp_summary <- x_wide %>%
  group_by(inhib_const_mult, Species) %>%
  summarise(hyst_1 = mean(abs(log10(up) - log10(down))),
            hyst_2 = get_hysteresis_range(up, down, a),
            nl_up = get_nonlinearity(a, log10(up)),
            nl_down = get_nonlinearity(a, log10(down))) %>%
  ungroup() %>%
  gather(key = resp_measure, value = value, 3:6) %>%
  mutate(
         var_type = ifelse(str_sub(Species, 1, 1)=="N", "Organism", "Substrate")
         )

6.3.4.1 All species nonlinearity

resp_summ1 <- resp_summary %>%
  filter(resp_measure %in% c("nl_down", "nl_up")) %>%
  group_by(inhib_const_mult, Species, var_type) %>%
  summarise(ave_nl = mean(value))

resp_summ1 %>%
  ggplot() +
  geom_point(aes(inhib_const_mult, ave_nl)) +
  facet_wrap( ~ Species)

6.3.4.2 Substrate nonlinearity

focus_on <- "Substrate"
ave_nl_substrate <- resp_summ1 %>%
  filter(var_type == focus_on) %>%
  group_by(inhib_const_mult) %>%
  summarise(ave_ave_nl = mean(ave_nl))
ave_nl_substrate %>%
  ggplot() +
  geom_point(aes(inhib_const_mult, ave_ave_nl)) +
  ggtitle(paste("Nonlinearity in", focus_on, "response."))

6.3.4.3 Organism nonlinearity

focus_on <- "Organism"
ave_nl_organism <- resp_summ1 %>%
  filter(var_type == focus_on) %>%
  group_by(inhib_const_mult) %>%
  summarise(ave_ave_nl = mean(ave_nl))
ave_nl_organism %>%
  ggplot() +
  geom_point(aes(inhib_const_mult, ave_ave_nl)) +
  ggtitle(paste("Nonlinearity in", focus_on, "response."))

6.3.4.4 System nonlinearity

focus_on <- "System"
ave_nl_system <- resp_summ1 %>%
  group_by(inhib_const_mult) %>%
  summarise(ave_ave_nl = mean(ave_nl))
ave_nl_system %>%
  ggplot() +
  geom_point(aes(inhib_const_mult, ave_ave_nl)) +
  ggtitle(paste("Nonlinearity in", focus_on, "response."))

ggsave("figures/system_nonlinearity.png", width = 3, height = 2)

6.3.5 Hysteresis

6.3.5.1 Substrate hysteresis

focus_on <- "Substrate"
resp_summ3 <- resp_summary %>%
  filter(resp_measure %in% c("hyst_1", "hyst_2"),
         var_type == focus_on) %>%
  group_by(inhib_const_mult, resp_measure, var_type) %>%
  summarise(ave_hyst = mean(value))

resp_summ3 %>%
  ggplot() +
  geom_point(aes(x = inhib_const_mult, y = ave_hyst)) +
  facet_wrap( ~ resp_measure, scales="free") +
  ggtitle(paste("Hysteresis in", focus_on, "response."))

6.3.5.2 Organism hysteresis

focus_on <- "Organism"
resp_summ3 <- resp_summary %>%
  filter(resp_measure %in% c("hyst_1", "hyst_2"),
         var_type == focus_on) %>%
  group_by(inhib_const_mult,  resp_measure, var_type) %>%
  summarise(ave_hyst = mean(value))

resp_summ3 %>%
  ggplot() +
  geom_point(aes(x = inhib_const_mult, y = ave_hyst)) +
  facet_wrap( ~ resp_measure, scales="free") +
  ggtitle(paste("Hysteresis in", focus_on, "response."))

6.3.5.3 System hysteresis (separate)

resp_summ4 <- resp_summary %>%
  filter(resp_measure %in% c("hyst_1", "hyst_2")) %>%
  group_by(resp_measure, inhib_const_mult) %>%
  summarise(mean_hyst = mean(value))

resp_summ4 %>%
  ggplot() +
  geom_point(aes(x = inhib_const_mult, y = mean_hyst)) +
  facet_wrap( ~ resp_measure) +
  ggtitle("Hysteresis in system response.")

6.3.5.4 System hysteresis (combined)

resp_summ5 <- resp_summ4 %>%
  group_by(inhib_const_mult) %>%
  summarise(mean_mean_hyst = mean(mean_hyst))

resp_summ5 %>%
  ggplot() +
  geom_point(aes(x = inhib_const_mult, y = mean_mean_hyst)) +
  xlab("Organismal tolerance") +
  ylab("Amount of hysteresis") +
  theme_light() +
  ylim(0,2)

ggsave("figures/inhibition_system_hysteresis.png", width = 2.5, height = 2)

7 Results: Temporal state switching

7.1 Switch anoxic-oxic-anoxic

## anoxic initial state
init_state <- init_state_favouring_anoxic

## simulation timing
simulation_duration <- 5000
simulation_sampling_interval <- 10
event_interval_t <- 100

## low-high-low oxgygen diffusivity 
lhl_log10a_series <- c(-8, -8, -2, -2, -2, -2, -2, -8, -8)



## run simulation
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic,
                          log10a_series = lhl_log10a_series,
                          sim_duration = simulation_duration,
                          sim_sample_interval = simulation_sampling_interval,
                          event_interval = event_interval_t)
saveRDS(sim_res, "data/lhl_with_defaults.RDS")
lhl_with_defaults <- readRDS("data/lhl_with_defaults.RDS")
plot_dynamics(lhl_with_defaults)

lhl_with_defaults$result %>%
  filter(between(time, 1000, 73000)) %>%
  ggplot() +
  geom_path(aes(x=a, y=log10(O)),
            arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")))

7.2 Switch oxic-anoxic-oxic

## anoxic initial state
init_state <- init_state_favouring_oxic

## simulation timing
simulation_duration <- 50000
simulation_sampling_interval <- 10
event_interval_t <- 100

## low-high-low oxgygen diffusivity 
hlh_log10a_series <- c(-2, -2, -10, -14, -18, -22, -2, -2, -2)



## run simulation
sim_res <- run_simulation(initial_state = init_state_favouring_oxic,
                          log10a_series = hlh_log10a_series,
                          sim_duration = simulation_duration,
                          sim_sample_interval = simulation_sampling_interval,
                          event_interval = event_interval_t)
saveRDS(sim_res, "data/hlh_with_defaults.RDS")
hlh_with_defaults <- readRDS("data/hlh_with_defaults.RDS")
plot_dynamics(hlh_with_defaults)

hlh_with_defaults$result %>%
  filter(between(time, 1000, 100000)) %>%
  ggplot() +
  geom_path(aes(x=a, y=log10(O)),
            arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")))

8 Appendices

8.1 runsteady function not finding correct values

init_state_anoxic_2 <- c(N_CB = 1e1,
           N_PB = 1e8,
           N_SB = 1e8,
           SO = 300,
           SR = 300,
           O = 1e-5,
           P = 1e1)


RS <- runsteady(y = init_state_anoxic_2,
                fun = bushplus_dynamic_model,
                parms = default_parameter_values,
                times = c(0, 1e5),
                stol=1)
RS$y
RS <- runsteady(y = init_state_favouring_oxic,
                fun = bushplus_dynamic_model,
                parms = default_parameter_values,
                times = c(0, 1e5),
                stol=1)
RS$y

8.2 Separatrix

Note that this section not working / updated since the code was made into functions

8.2.1 Effect of N_CB and oxygen diffusability

noise_sigma <- 0.0

get_final_states_N_CB <- function(x) { 
  
  #give N_CB, N_SB, N_PB and a_O values and put these into state and parameters
  state["N_CB"] <- x["N_CB"]
  parameters["a_O"] <- x["a_O"]
  
  ## update forcing function
  a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0,
                                             log10(parameters["a_O"]),
                                             max(times),
                                             log10(parameters["a_O"])))
  log10a_forcing2 <- approxfun(x = a_forcing1[,1],
                             y = a_forcing1[,2],
                             method = "linear",
                             rule = 2)
  
  ## assign the changed a_forcing2 into the global environment, from where the model gets a_forcing2
  assign("log10a_forcing2", log10a_forcing2, envir = .GlobalEnv)

  as.data.frame(ode(y = state,
                    times = times,
                    func = micro_mod_var,
                    parms = parameters,
                    method = "radau",
                    events = list(func=noise_events, time=times)))[2,-1] 
  #gives you values for the last time step
}

times <- seq(0, 50000, by = 50000)

state <- c(N_CB = NA, #this one is varied
           N_PB = 1e8,
           N_SB = 1e8,
           SO = 200,
           SR = 200,
           O = 10,
           P = 9.5)

#different combinations of N_CB and oxygen diffusability
grid_num <- 50
expt <- expand.grid(N_CB = 10^seq(0, 8, length=grid_num),
                    a_O = 10^seq(-3.5, -2, length=grid_num))


result <- apply(expt, 1, function(x) get_final_states_N_CB(x)) %>%
                          tibble() %>%
                          unnest() %>%
                          mutate(initial_N_CB=expt$N_CB,
                                 a_O=expt$a_O,
                                 highlow=ifelse(N_CB>2e8, 'high', 'low'))

result <- result %>%
  mutate(N_CB = ifelse(N_CB < minimum_abundance, minimum_abundance, N_CB))



result$highlow <- as.factor(result$highlow)
saveRDS(result, "data/CB_oxdiff_hyst.RDS")
result <- readRDS("data/CB_oxdiff_hyst.RDS")
#plot N_CB / P_CB against a_0
p <- ggplot(result) + 
  geom_line(aes(x=log10(a_O), y=log10(N_CB), color=highlow), size=1.25) +
  ggtitle("Stable states of the system")+
  xlab('log[oxygen diffusivity]') +
  ylab('Initial cyanobacteria abundance')
p <- p + geom_point(data=result, aes(x=log10(a_O), y=log10(initial_N_CB), color=highlow)) + labs(color = "Final CB concentration")
plot(p)
#ggsave("initial conditions_1.png")
#show hysterisis effects for all the parameters
#plot oxygen against diffusability
ggplot(result) + 
  geom_point(aes(x=log10(a_O), y=(log10(O))))+
  xlab('log[oxygen diffusivity]') +
  ylab('log[oxygen concentration]')+
  geom_segment(aes(x=log10(result$a_O[839]), y=log10(result$O[839]), xend=log10(result$a_O[838]), yend=log10(result$O[838])), arrow=arrow(), colour="blue", size=0.5)+
  geom_segment(aes(x=log10(result$a_O[1761]), y=log10(result$O[1761]), xend=log10(result$a_O[1762]), yend=log10(result$O[1762])), arrow=arrow(), colour="blue", size=0.5)+
  geom_segment(aes(x=-4.25, y=-0.45, xend=-3.75, yend=0.1), arrow=arrow(), colour="blue", size=0.5)+
  geom_curve(aes(x = -3.75, y = 1.9, xend = -4.25, yend = 1.5), arrow=arrow(), colour='blue', size = 0.5, curvature = -0.2)
#ggsave("hysterisis_CB_1.png")
ggplot(result) + 
  geom_point(aes(x=log10(a_O), y=(log10(P))))+
  ggtitle("Hysterisis effect of phoshate concentration")+
  xlab('log[oxygen diffusivity]') +
  ylab('log[phosphate concentration]')
#ggsave("hysterisis_CB_2.png")
ggplot(result) + 
  geom_point(aes(x=log10(a_O), y=(log10(SO))))+
  ggtitle("Hysterisis effect of oxidized sulfur concentration")+
  xlab('log[oxygen diffusivity]') +
  ylab('log[oxydized sulfur concentration]')
#ggsave("hysterisis_CB_3.png")
ggplot(result) + 
  geom_point(aes(x=log10(a_O), y=(log10(SR))))+
  ggtitle("Hysterisis effect of reduced sulfur concentration")+
  xlab('log[oxygen diffusivity]') +
  ylab('log[reduced sulfur] concentration]')
#ggsave("hysterisis_CB_4.png")

These graphs show an at least qualitative match to the result in the paper.

8.2.2 Effect of PB/SB depending on oxygen diffusability

noise_sigma <- 0.00
get_final_states_N_PB <- function(x) { 
  #give in N_CB, N_SB, N_PB and a_O values and put these into state and parameters
  state["N_PB"] <- x["var"]
  state["N_SB"] <- x["var"]
  parameters["a_O"] <- x["a_O"]
  
  ## update forcing function
  a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0,
                                             log10(parameters["a_O"]),
                                             max(times),
                                             log10(parameters["a_O"])))
  log10a_forcing2 <- approxfun(x = a_forcing1[,1],
                             y = a_forcing1[,2],
                             method = "linear",
                             rule = 2)
  ## assign the changed a_forcing2 into the global environment, from where the model gets a_forcing2
  assign("log10a_forcing2", log10a_forcing2, envir = .GlobalEnv)
  
  as.data.frame(ode(y = state,
                         times = times,
                         func = micro_mod_var,
                         parms = parameters,
                         method = "radau", events = list(func=noise_events, time=times)))[2,-1] 
  #gives you values for the last time step
}

times <- seq(0, 50000, by = 5000)
state <- c(N_CB = 5, 
           N_PB = NA,#this one is varied
           N_SB = NA,
           SO = 200,
           SR = 200,
           O = 10,
           P = 9.5)

#different combinations of N_PB and oxygen diffusivity
grid_num <- 50
expt <- expand.grid(var=10^seq(0, 7, length=grid_num),
                    a_O = 10^seq(-3, -2, length=grid_num))


result2 <- apply(expt, 1, function(x) get_final_states_N_PB(x)) %>%
                          tibble() %>%
                          unnest() %>%
                          mutate(initial_N_PB=expt$var,
                                 initial_N_SB=expt$var,
                                 a_O=expt$a_O,
                                 highlow=ifelse(N_PB>1e8, 'high', 'low')) 

result2 <- result2 %>%
  mutate(N_PB = ifelse(N_PB < minimum_abundance, minimum_abundance, N_PB))

result2$highlow <- as.factor(result2$highlow)
saveRDS(result2, "data/PBSB_oxdiff_hyst.RDS")
result2 <- readRDS("data/PBSB_oxdiff_hyst.RDS")
#plot N_PB against a_0
p <- ggplot(result2) + 
  geom_line(aes(x=log10(a_O), y=log10(N_PB), color=highlow), size=1.25) +
  ggtitle("Stable states of the system")+
  xlab('log[oxygen diffusivity]') +
  ylab('Initial PS and SR bacteria concentration')
p <- p + geom_point(data=result2,
                    aes(x=log10(a_O), y=log10(initial_N_PB), color=highlow))
p <- p + labs(color = "Final PB and SB concentration")
plot(p)
#ggsave("initial conditions_2.png")

8.3 EWS analyses

Evaluation of the switch: Get the functions:

source('C:/Users/probs/Desktop/Sem 6/Internship/functions_analysis_EWS_2.R')

Get state shifts and calculate parameters:

x3 <- 30 #size of wished rolling window (percentage btw 0 and 100)
x5 <- 50 #size of wished rolling window (percentage btw 0 and 100)

#first switch
name_1='anoxic_oxic 1'
trans1 <- data.frame(out_t %>% filter(time>87500 & time<92500)) #for analysis
trans1_f <- data.frame(out_f %>% filter(time>87500 & time<92500)) #for plotting
n_13 <- round(length(trans1$time)*x3/100) #size of rolling window
t_13 <- length(trans1$time) - n_13 + 1 #number of rolling windows
n_15 <- round(length(trans1$time)*x5/100) #size of rolling window
t_15 <- length(trans1$time) - n_15 + 1 #number of rolling windows
out_df_1_N_CB <- data.frame(trans1$N_CB)
colnames(out_df_1_N_CB) <- c("N_CB")
out_df_1_N_PB <- data.frame(trans1$N_PB)
colnames(out_df_1_N_PB) <- c("N_PB")
out_df_1_N_SB <- data.frame(trans1$N_SB)
colnames(out_df_1_N_SB) <- c("N_SB")

ggplot(aes(x=time, y=trans_quantity, col=species), data=trans1_f) +
  ggtitle("Anoxic_oxic 1")+
  xlab('Time') +
  ylab('Trans_quantity') +
  facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
  geom_line()+
  labs(color = "Organism/Substrate")
ggsave("anoxic_oxic 1.png")

#second switch
name_2='oxic_anoxic 1'
trans2 <- data.frame(out_t %>% filter(time>335000 & time<360000))
trans2_f <- data.frame(out_f %>% filter(time>335000 & time<360000))
n_23 <- round(length(trans2$time)*x3/100) 
t_23 <- length(trans2$time) - n_23 + 1
n_25 <- round(length(trans2$time)*x5/100) 
t_25 <- length(trans2$time) - n_25 + 1
out_df_2_N_CB <- data.frame(trans2$N_CB)
colnames(out_df_2_N_CB) <- c("N_CB")
out_df_2_N_PB <- data.frame(trans2$N_PB)
colnames(out_df_2_N_PB) <- c("N_PB")
out_df_2_N_SB <- data.frame(trans2$N_SB)
colnames(out_df_2_N_SB) <- c("N_SB")

ggplot(aes(x=time, y=trans_quantity, col=species), data=trans2_f) +
  ggtitle("Oxic_anoxic 1")+
  xlab('Time') +
  ylab('Trans_quantity') +
  facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
  geom_line()+
  labs(color = "Organism/Substrate")
ggsave("oxic_anoxic 1.png")

#third switch
name_3='oxic_anoxic 2'
trans3 <- data.frame(out2_t %>% filter(time>165000 & time<185000))
trans3_f <- data.frame(out2_f %>% filter(time>165000 & time<185000))
n_33 <- round(length(trans3$time)*x3/100) 
t_33 <- length(trans3$time) - n_33 + 1
n_35 <- round(length(trans3$time)*x5/100) 
t_35 <- length(trans3$time) - n_35 + 1
out2_df_1_N_CB <- data.frame(trans3$N_CB)
colnames(out2_df_1_N_CB) <- c("N_CB")
out2_df_1_N_SB <- data.frame(trans3$N_SB)
colnames(out2_df_1_N_SB) <- c("N_SB")
out2_df_1_N_PB <- data.frame(trans3$N_PB)
colnames(out2_df_1_N_PB) <- c("N_PB")

ggplot(aes(x=time, y=trans_quantity, col=species), data=trans3_f) +
  ggtitle("Oxic_anoxic 2")+
  xlab('Time') +
  ylab('Trans_quantity') +
  facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
  geom_line()+
  labs(color = "Organism/Substrate")
ggsave("oxic_anoxic 2.png")

#fourth switch
name_4='anoxic_oxic 2'
trans4 <- data.frame(out2_t %>% filter(time>340000 & time<348500))
trans4_f <- data.frame(out2_f %>% filter(time>340000 & time<348500))
n_43 <- round(length(trans4$time)*x3/100) 
t_43 <- length(trans4$time) - n_43 + 1
n_45 <- round(length(trans4$time)*x5/100) 
t_45 <- length(trans4$time) - n_45 + 1
out2_df_2_N_CB <- data.frame(trans4$N_CB)
colnames(out2_df_2_N_CB) <- c("N_CB")
out2_df_2_N_PB <- data.frame(trans4$N_PB)
colnames(out2_df_2_N_PB) <- c("N_PB")
out2_df_2_N_SB <- data.frame(trans4$N_SB)
colnames(out2_df_2_N_SB) <- c("N_SB")

ggplot(aes(x=time, y=trans_quantity, col=species), data=trans4_f) +
  ggtitle("Anoxic_oxic 2")+
  xlab('Time') +
  ylab('Trans_quantity') +
  facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
  geom_line()+
  labs(color = "Organism/Substrate")
ggsave("anoxic_oxic 2.png")

Rolling window analysis

#analysis for 1st
rollingwindow_analysis(trans1, 2, n_15, t_15, name_1)
rollingwindow_analysis(trans1, 3, n_15, t_15, name_1)
rollingwindow_analysis(trans1, 4, n_15, t_15, name_1)
  
#analysis for 2nd
rollingwindow_analysis(trans2, 2, n_25, t_25, name_2)
rollingwindow_analysis(trans2, 3, n_25, t_25, name_2)
rollingwindow_analysis(trans2, 4, n_25, t_25, name_2)

#analysis for 3rd
rollingwindow_analysis(trans3, 2, n_35, t_35, name_3)
rollingwindow_analysis(trans3, 3, n_35, t_35, name_3)
rollingwindow_analysis(trans3, 4, n_35, t_35, name_3)

#analysis for 4th
rollingwindow_analysis(trans4, 2, n_45, t_45, name_4)
rollingwindow_analysis(trans4, 3, n_45, t_45, name_4)
rollingwindow_analysis(trans4, 4, n_45, t_45, name_4)

Metric based indicators from the earlywarnings package

#analysis for 1st
earlywarnings_plot(out_df_1_N_CB, name_1, x3)
earlywarnings_plot(out_df_1_N_PB, name_1, x3)
earlywarnings_plot(out_df_1_N_SB, name_1, x3)

#analysis for 2nd
earlywarnings_plot(out_df_2_N_CB, name_2, x3)
earlywarnings_plot(out_df_2_N_PB, name_2, x3)
earlywarnings_plot(out_df_2_N_SB, name_2, x3)

#analysis for 3rd
earlywarnings_plot(out2_df_1_N_CB, name_3, x3)
earlywarnings_plot(out2_df_1_N_PB, name_3, x3)
earlywarnings_plot(out2_df_1_N_SB, name_3, x3)

#analysis for 4th
earlywarnings_plot(out2_df_2_N_CB, name_4, x3)
earlywarnings_plot(out2_df_2_N_PB, name_4, x3)
earlywarnings_plot(out2_df_2_N_SB, name_4, x3)

Conditional Heteroskedasticity

#analysis for 1st
ch <- condhet_analysis(out_df_1_N_CB, name_1)
ch <- condhet_analysis(out_df_1_N_PB, name_1)
ch <- condhet_analysis(out_df_1_N_SB, name_1)

#analysis for 2nd
ch2 <- condhet_analysis(out_df_2_N_CB, name_2)
ch2 <- condhet_analysis(out_df_2_N_PB, name_2)
ch2 <- condhet_analysis(out_df_2_N_SB, name_2)

#analysis for 3rd
ch3 <- condhet_analysis(out2_df_1_N_CB, name_3)
ch3 <- condhet_analysis(out2_df_1_N_PB, name_3)
ch3 <- condhet_analysis(out2_df_1_N_SB, name_3)

#analysis for 4th
ch4 <- condhet_analysis(out2_df_2_N_CB,name_4)
ch4 <- condhet_analysis(out2_df_2_N_PB, name_4)
ch4 <- condhet_analysis(out2_df_2_N_SB, name_4)

Detrending fluctuation analysis

#analysis for 1st
dfa_analysis(trans1, 2, n_15, t_15, name_1)
dfa_analysis(trans1, 3, n_15, t_15, name_1)
dfa_analysis(trans1, 4, n_15, t_15, name_1)

#analysis for 2nd
dfa_analysis(trans2, 2, n_25, t_25, name_2)
dfa_analysis(trans2, 3, n_25, t_25, name_2)
dfa_analysis(trans2, 4, n_25, t_25, name_2)

#analysis for 3rd
dfa_analysis(trans3, 2, n_35, t_35, name_3)
dfa_analysis(trans3, 3, n_35, t_35, name_3)
dfa_analysis(trans3, 4, n_35, t_35, name_3)

#analysis for 4th
dfa_analysis(trans4, 2, n_45, t_45, name_4)
dfa_analysis(trans4, 3, n_45, t_45, name_4)
dfa_analysis(trans4, 4, n_45, t_45, name_4)

Drift Diffusion Jump Nonparametrics Early Warning Signals:

#shift 1
ddj11<- ddj_model(out_df_1_N_CB)
plot_ddj_1(out_df_1_N_CB, ddj11, name_1)
plot_ddj_2(out_df_1_N_CB, ddj11, name_1)

ddj12<- ddj_model(out_df_1_N_PB)
plot_ddj_1(out_df_1_N_PB, ddj12, name_1)
plot_ddj_2(out_df_1_N_PB, ddj12, name_1)

ddj13<- ddj_model(out_df_1_N_SB)
plot_ddj_1(out_df_1_N_SB, ddj13, name_1)
plot_ddj_2(out_df_1_N_SB, ddj13, name_1)


#shift 2
ddj21<- ddj_model(out_df_2_N_CB)
plot_ddj_1(out_df_2_N_CB, ddj21, name_2)
plot_ddj_2(out_df_2_N_CB, ddj21, name_2)

ddj22<- ddj_model(out_df_2_N_PB)
plot_ddj_1(out_df_2_N_PB, ddj22, name_2)
plot_ddj_2(out_df_2_N_PB, ddj22, name_2)

ddj23<- ddj_model(out_df_2_N_SB)
plot_ddj_1(out_df_2_N_SB, ddj23, name_2)
plot_ddj_2(out_df_2_N_SB, ddj23, name_2)


#shift 3
ddj31<- ddj_model(out2_df_1_N_CB)
plot_ddj_1(out2_df_1_N_CB, ddj31, name_3)
plot_ddj_2(out2_df_1_N_CB, ddj31, name_3)

ddj32<- ddj_model(out2_df_1_N_PB)
plot_ddj_1(out2_df_1_N_PB, ddj32, name_3)
plot_ddj_2(out2_df_1_N_PB, ddj32, name_3)

ddj33<- ddj_model(out2_df_1_N_SB)
plot_ddj_1(out2_df_1_N_SB, ddj33, name_3)
plot_ddj_2(out2_df_1_N_SB, ddj33, name_3)


#shift 4
ddj41<- ddj_model(out2_df_2_N_CB)
plot_ddj_1(out2_df_2_N_CB, ddj41, name_4)
plot_ddj_2(out2_df_2_N_CB, ddj41, name_4)

ddj42<- ddj_model(out2_df_2_N_PB)
plot_ddj_1(out2_df_2_N_PB, ddj42, name_4)
plot_ddj_2(out2_df_2_N_PB, ddj42, name_4)

ddj43<- ddj_model(out2_df_2_N_SB)
plot_ddj_1(out2_df_2_N_SB, ddj43, name_4)
plot_ddj_2(out2_df_2_N_SB, ddj43, name_4)